Unified transient creep constitutive model based on the crack evolution of micritic bioclastic limestone

The surrounding rock at the exit of the No. 1 drainage tunnel of the Artashi Water Conservancy Project is micritic bioclastic limestone with 55% bioclastic material. This rock underwent unpredictable large and time-dependent deformation during excavation. To date, the mechanical behaviour of this kind of rock has rarely been studied. In this study, traditional triaxial compression tests and multilevel creep tests were conducted on micritic bioclastic limestone, and the results clarified the instantaneous and time-dependent mechanical properties of the rock. Considering that the essence of rock failure is crack growth, the crack strain evolution properties were revealed in rock triaxial compression tests and multilevel creep tests. Based on triaxial compression tests, the evolution of axial cracks with increasing deviatoric stress ratio Rd (ratio of deviatoric stress to peak deviatoric stress) was observed, and an axial crack closure element and new crack growth element were proposed. To simulate the creep behaviour of a rock specimen, the relationship of the rock creep crack strain rate with Rd was studied. A creep crack element was created, and the creep crack strain evolution equation was obtained, which closely fit the experimental data. Combining the 4 element types (elastic element, crack closure element, crack growth element, and creep crack element), a unified transient creep constitutive model (Mo’s model) was proposed, which represented both the transient and time-dependent mechanical properties of the micritic bioclastic limestone.


Introduction
The Artashi Water Conservancy Project in Xinjiang Province, China, is the controlling project on the Yarkant River, undertaking flood control, irrigation diversion and power generation. The surrounding rock at the exit of No. 1 drainage tunnel of the Artashi Water Conservancy Project has experienced unexpected large and continuous deformation, leading to safety problems of the tunnel and slope. The deformation of the surrounding rock and slope have increased with time, resulting in fracture of the slope concrete and yield failure of the steel a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 arch of the tunnel (Fig 1(A)). The relative displacement of the outer surface of the arch top recorded by a 4-point monitor (M4-01) was 93.77 mm after 1332 days. The slope and the tunnel exit were reinforced with anchor cable and concrete anchor supports. Then monitoring and on-site observations revealed that slope deformation continued and the support structure underwent cracking failure, as shown in Fig 1(B). Therefore, it is necessary to investigate the composition and mechanical properties of the rock to reveal the mechanisms of rock deformation.
The surrounding rock is black and has a relatively dense massive structure, and thin section analysis results show that the rock has a biological skeleton structure (Fig 2). Numerous incipient cracks can also be seen in the rock. The composition of this rock is 45% interstitial material and 55% bioclastic material, occasionally with terrigenous clastic material. The widely distributed interstitial material is mainly a micrite substrate composed of micrite calcite. Bioclasts are mostly round, oval, tubular, curved and irregular in shape, and the grains vary in size, with diameters mostly less than 2 mm. The particle diameter size of the terrigenous clastic material is generally smaller than 0.05 mm, mainly quartz, and sporadically distributed. Accordingly, this rock is named micritic bioclastic limestone.
To the authors' knowledge, there are very few reports of in-depth studies of this kind of rock. This type of rock was previously encountered in engineering, mainly near faults or weak structural surfaces, and rarely with bioclastic limestone as the tunnel surrounding rock. Studies on limestone containing bioclastic material have been conducted by many scholars. EL-Sorogy AS et al. [1] performed a study on the composition and diagenesis of reefal limestone in central Saudi Arabia. Rožič B et al. [2] studied the composition of Middle Jurassic limestone megabreccia from the Slovenian Basin and inferred the evolution of the platform margin. V. Vajdova et al. [3,4] used microscopy to study microstructural damage and pore changes in porous carbonates containing bioclasts after triaxial compression experiments at different confining pressures; the results revealed that the microscopic damage mechanism was pore collapse and crack coalescence. Ji. Y et al. [5] used X-ray microtomography 3D imaging to study the pore changes in Indiana limestone containing bioclastic material after compression experiments and observed the pore reduction characteristics. The mechanical properties of reef limestone, which also contained bioclasts, have been progressively studied in recent years. Changqi Z et al. [6] investigated the microstructures and uniaxial compressive strength of beach calcarenites in the South China Sea. Q S. Meng et al. [7] conducted Hopkinson pressure bar tests to study the dynamic failure properties of coral reef limestones. The crack initiation and damage evolution of micritized framework reef limestone from the South China Sea were studied by triaxial compression tests [8]. However, research on the time-dependent mechanical properties of micritic bioclastic limestone has rarely been reported. As mentioned by H. Liu et al. [8], the mechanical properties of reef limestone depend on their composition, structure and diagenesis. The complex composition and specific structure of this micritic bioclastic limestone make it very unique, and thus, it is valuable to study the mechanical properties of this rock.
The phenomenon of time-dependent deformation of the No. 1 drainage tunnel surrounding rock and slope is a typical creep-type deformation [9][10][11]. Many scholars have conducted numerous creep tests, such as single-stage, multilevel creep tests and loading cycle creep tests, on different rocks to investigate their creep behaviour [9,[11][12][13][14][15]. To characterize the creep evolution of rocks, empirical models, component models, and other novel models have been proposed. An empirical model can represent the creep properties of a particular rock accurately with relatively few parameters. The drawback of an empirical model is that the parameters always lack physical meaning. Empirical models are only suitable for describing creep test data at specific stresses and is difficult to apply them to make engineering predictions for complex stress conditions [16][17][18][19]. For a component model, a spring, dashpot and slider are used to represent the elastic, viscous and plastic mechanical response of rock, respectively. With the combination of elements, this model can represent the primary and secondary creep stages, such as the models of Maxwell, Kelvin, Burgers, and Nishihara [20][21][22][23][24]. Fractional order has been introduced into the viscous element to represent the tertiary creep stage [25][26][27][28][29]. A component model is relatively easy to understand, and the components and parameters have physical meaning. However, a component model is phenomenologically proposed based on the mechanical response of rock and not based on the mechanism of creep [11,30].
To study the creep properties of rocks in depth and to develop constitutive models, important work has been performed by a number of scholars in recent years. A novel damage-based creep model has been proposed based on the time-hardening and damage theory of rock [30]. Zhang Xing et al. [31] established the damage evolution equation based on the variation in the deformation modulus of rock in multistage creep experiments, and then the creep equation was established. Based on the growth model of a single crack, scholars have used the micro/ macro approach to establish the relationship between crack extension and macrostrain and proposed the brittle creep model [32]. Yanlin Zhao et al. [33] Based on wing crack propagation theory, crack propagation properties were investigated, and an equivalent Burgers model was proposed to characterise the normalised crack propagation length. With the development of particle flow codes, numerical methods have been adopted to investigate the creep properties of rocks. Hua Li et al. [34] proposed a 3D microcreep model to simulate the creep properties of rock salt based on PFC3D, which agreed well with indoor tests. Kefeng Zhou et al. [35] investigated the effect of shear parameters on the creep properties of rocks containing joints based on PFC numerical software. The microscopic mechanism of rock creep and damage has been generally accepted as the evolution of microcracks. The crack strain proposed by Martin [36] ignored the complex and random distribution of initial cracks and thus avoided the difficulties caused by uncertainty [37]. In the present study, micritic bioclastic limestone was obtained from the Artashi Water Conservancy Project, and conventional triaxial and multilevel loading creep tests were conducted. The crack strain characteristics of rock specimens in triaxial compression and creep tests were revealed, and the crack evolution at all stages of the creep process was investigated. A unified transient creep constitutive model based on crack evolution was proposed.

Rock experiment preparation and procedures
The rock used to make the specimens was taken from the Artashi Water Conservancy Project in Xinjiang, China. A whole rock block, approximately 900 kg, was obtained from the surrounding rock of the large-deformation tunnel to reduce variability among the rock specimens. All rock samples were drilled in the same direction and carefully processed into standard cylindrical rock specimens with dimensions of F50 mm×100 mm (Fig 3) and a height-to-diameter ratio of 2.0. The average density of the rock specimens is 2.77 g/cm 3 .
The confining pressures adopted in the conventional triaxial compression testing were 1, 3, 5 and 10 MPa, with 5 specimens tested under each condition. The loading equipment in these tests was the MTS815 apparatus at Sichuan University. The specifications of the loading equipment are shown in Table 1 [38].
The conventional triaxial compression tests of the micritic bioclastic limestone were performed as follows: (1) The confining pressure was applied to a predetermined value with a loading rate of 6 MPa/min.
(2) An axial load was applied at a loading rate of 20 kN/min.
(3) The axial load control method was changed to lateral deformation rate control with 0.02 mm/min at the end of the elastic stage of the rock specimen.
(4) After the peak stress, the axial load was applied at a lateral deformation rate of 0.04 mm/ min until the residual deformation stage of the rock was reached.
The multilevel loading creep tests on the micritic bioclastic limestone were carried out in the following steps: (1) The confining pressure and axial load were applied in the same way as in the conventional triaxial compression test.
(2) After the axial load was applied to a predetermined value, the deviatoric stress was kept constant for a predetermined time (2 h).
(3) After the creep period ended, the load was unloaded at a rate of 20 kN/min until the deviatoric stress reached 0.25 MPa, which was held for 0.2 h.
(4) The next deviatoric stress level was applied, and steps (2) and (3) were repeated until the specimen failed.
The confining pressure remained constant for each level of the multilevel loading creep tests. Based on the stress thresholds of the triaxial compression test, as a load level reference, 4

Laboratory test of micritic bioclastic limestone
To study the time-dependent deformation and crack evolution characteristics of the tunnels, conventional triaxial and multilevel loading creep test experiments on the rock were carried out. The results of the triaxial compression tests are shown in Based on the results of the triaxial compression tests, the stress thresholds of micritic bioclastic limestone are determined in Section 3. Considering the rock strain response during the test, the adjusted deviatoric stress of multilevel loading creep tests are shown in Table 2.
The multilevel loading rock creep tests were carried out by MTS815, and the results are shown in Fig 5. As the load reaches the predetermined value, the axial strain of the rock evolves with time, exhibiting significant creep characteristics. The axial strain and the slope of the curves increase with the loading level.

Crack strain calculation and stress threshold determination
Axial and lateral strains arise when rock specimens are subjected to forces, including elastic and crack strains. The elastic strain of loaded specimens in a conventional triaxial compression test can be calculated by Hooke's law and can be denoted by ε 1e and ε 3e in the axial and lateral directions: where ε 1e and ε 3e denote the axial and lateral elastic strains, respectively, E is the elastic modulus, μ e is the elastic Poisson's ratio, σ 1 denotes the axial stress, and σ 3 is the confining pressure. From Eq (1), the axial strain increment independent of Poisson's ratio can be expressed as: According to Reference [39], crack strain refers to the strain resulting from the closure, propagation and coalescence of cracks. The crack strain of the rock specimens includes the axial crack strain ε 1c and the lateral crack strain ε 3c . The axial and lateral strains obtained from triaxial compression tests consist of elastic and crack strains, which are presented in Eq (3): ( where ε 1 denotes the axial strain and ε 3 denotes the lateral strain. Therefore, the rock crack strain can be obtained from Eq (4): The rock crack strain in Eq (4) is sensitive to the elastic parameters (E and μ e ). Thus, E and μ e of the rock specimen must be obtained accurately. The elastic modulus E of rock specimens under triaxial compression can be obtained by the slope of the linear portion of the axial stress-strain curve [40,41]. A method called the moving point regression technique, which uses a sliding window approach to move the axial strain-stress dataset to fit a straight line over an interval, has been proposed to obtain the elastic modulus of rocks [42]. According to Reference [42], a relatively smooth tangent modulus curve can be obtained when 5% of the total number of data points is within a certain interval (Fig 6). Eq (1) indicates that the lateral elastic strain depends on the elastic Poisson's ratio. To accurately determine the elastic Poisson's ratio of rock, the moving point regression technique is adopted.
Reference [44] defined the failure process of brittle rock in five stages by four stress thresholds: crack closure stress σ cc , crack initiation stress σ ci , crack damage σ cd stress and peak stress σ p . Reference [45] reported many methods for determining the stress thresholds of rock based on previous studies. However, the deviatoric stress-volumetric strain curve shown in Fig 4 of micritic bioclastic limestone under triaxial compression has no typical reversal point, so the method proposed in Reference [36,45] cannot be used in determining σ cd . Therefore, an optimized method to determine the stress thresholds of rocks is proposed as follows.  (1) and (4).
6. Plot the deviatoric stress against the lateral and volumetric crack strains.
7. Determine σ ci by the point at the end of the vertical section of the deviatoric stress-lateral crack strain curve.
8. Determine σ cd as the point at the end of the vertical section of the deviatoric stress-axial crack strain curve. In step (3), σ cc is determined by the axial crack strain instead of the volumetric crack strain, thus avoiding the summation of errors and improving the accuracy of the results. In step (8), the axial crack strain curve is used instead of the volumetric strain curve, overcoming the problem caused by the absence of a typical reversal arc and reversal point of the volumetric strain curve. The stress thresholds obtained by Martin's method [36] and the present method are listed in Table 3. The stress thresholds (σ cc and σ ci ) determined by the two methods are close, whereas σ cd obtained by the present method is more reasonable than Martin's method.

Crack evolution characteristics of triaxial compression tests and creep tests
The deviatoric stress-crack strain curves of the micritic bioclastic limestone at confining pressures of 1 MPa, 3 MPa, 5 MPa and 10 MPa are represented in Fig 7. Because of the bioclastic The elastic modulus of each rock specimen during each loading of the multilevel creep tests was obtained by the moving point regression technique and is plotted in Fig 8. The curves clearly show that the elastic modulus of the rock tends to increase and then decrease with the loading level. When the rock undergoes the first few stages of loading and creep, the rock microcracks and pores close, the rock weak structures are compressed, and the elastic modulus of the rock increases. As the loading and creep stages increase, new cracks accumulate, and the elastic modulus of the rock gradually decreases. According to the elastic modulus at each stage of loading, the axial crack strain curve of the rock loading and creep process can be obtained. As shown in Fig 9, the axial crack strain can be divided into a loading stage and creep stage. The crack evolution characteristics in the loading stage are similar to those of the triaxial compression test. The cracks in the rock gradually close as the deviatoric stress increases, and then the rock enters the elastic deformation state. As the axial loading increases, the axial crack strain deviates from a straight line, indicating that the deviatoric stress exceeds σ cd and that axial cracks form in the rock, for example, see Fig  9(B). When the load holds at a constant deviatoric stress, the rock enters a creep state, in which the rock crack strain increases with holding time.
Taking level 1 of the 1 MPa confining pressure as an example, the percentage of creep crack strain is 64.8% of the initial crack closure strain, indicating that the creep crack strain at each creep level contributes negligibly to the total crack strain. Therefore, a constitutive model of rock that only considers transient rock deformation and ignores rock creep is inadequate for predicting the deformation and stability of underground caverns.

Crack strain element and elastic-crack constitutive model
To improve the constitutive model of rock, the transient crack evolution (loading) and the time-dependent crack evolution (creep) are unified based on crack strain. The following basic assumptions are introduced to study the time-dependent deformation characteristics of the micritic bioclastic limestone and establish the constitutive model.
(1) As proposed by Martin [36], the rock strain consists of elastic strain and crack strain.
(2) The elastic modulus of the rock matrix remains constant during a loading-creepunloading process.
(3) The rock strain caused by loading is assumed to be independent of time and related only to the deviatoric stress level.
(4) As the elastic strain of the rock remains constant during creep, the creep strain is considered to be the time-dependent crack strain of the rock.
In triaxial compression tests, brittle rock failure is always related to crack evolution, according to a reference [44]. This means that the crack strain can characterize the mechanical response of the rock. Therefore, the characteristics of the crack strain evolution can be studied to establish a relationship between them, which can be used to predict the rock mechanical properties. Based on the characteristic stresses from the triaxial compression test, R cc d is taken as 0.23, and R ci d is taken as 0.36. As the initial crack varies to some extent from specimen to specimen, the initial crack closure strain ε cc 1c is normalized to the initial crack closure strain ratio Rc (Eq (5)) to obtain the initial crack closure law.
The axial initial crack closure strain ratios at R c for confining pressures of 1, 3, 5 and 10 MPa are shown in Fig 10(A). It can be found that R c shows the same evolutionary pattern for different confining pressures. The crack evolution at the initial crack closure stage can be accurately expressed by fitting Eq (6).
where k 0 is named the initial crack closure factor, which can be fitted for any test sample, reflecting the crack closure characteristics.
The crack strain appears in the axial direction of the rock when the load increases to R cd d , marking the beginning of the unstable crack growth stage. To study the axial crack evolution in this rock without considering the variability of the specimens, the axial crack strain is normalized by the ratio of axial crack strain ε 1c to peak axial crack strain ε p 1c , or: The curves for R c versus R d are plotted in Fig 10(B). R c begins to develop when the rock bearing state exceeds R cd d and reaches 1, while R d equals 1. The evolutionary characteristics of the curves are the same for different confining pressures. The new axial crack strain evolution of the rock can be expressed by fitting Eq (8).
where k 1 denotes the axial new crack growth factor and C denotes the crack evolution factor. To characterize the rock crack strain behaviour under compression, rock crack evolution elements are proposed, as shown in Fig 11. The element used to represent the transient initial crack closure of the rock is shown in Fig 11(A). Based on the crack evolution of the initial crack closure stage, the equation for the initial crack closure element can be expressed as: From Fig 11(C), when R cc d � R d � R cd d , no crack grows in the axial direction of the rock specimen. Thus, Eq (9) can also be used between R cc d and R cd d . Similarly, the new crack growth element is shown in Fig 12(B), and the equation can be written as: The role of the axial crack strain elements in the rock compression process can be seen according to Fig 11(C). The initial crack strain is represented by the initial crack closure element at the beginning of the rock loading until R cc d . While the rock bearing state is between R cc d and R cd d , no new crack strain is generated in the axial direction. Once the rock bearing state exceeds R cd d , the new axial crack element comes into play, and the axial crack strain develops rapidly until Rd reaches 1.
The axial crack strain results from the triaxial compression tests can be represented by crack evolution elements. The crack strains ε cc 1c and ε p 1c are listed in Table 4. The LM (Levenberg-Marquardt) algorithm was adopted to fit the axial crack strain data and obtain the element parameters (C, k 0 , k 1 ). By fitting, the proposed initial crack closure and new crack growth elements can well characterize the evolution of the crack strain in the rock (Fig 12).
According to Eq (3), the total strain consists of elastic strain and crack strain. The elastic element is used to represent the elastic behaviour of the rock. Combined with the crack closure and new crack growth elements, a schematic representation of the elastic-crack constitutive model can be drawn as shown in Fig 13: The elastic element follows Hooke's law and can be expressed as Eq (11).
Accordingly, the elastic-crack constitutive model can be expressed as:

Creep crack strain element
The creep crack strain of micritic bioclastic limestone is calculated by the crack strain equation in Eq (4). Axial creep crack strain-time curves are shown in Fig 14. The creep crack increment strain of the first level is larger than that of the later level. In addition, the curves show that there are noteworthy increases in crack strain between the first and second creep levels, indicating that the initial cracks are significantly adjusted in the first creep level. The typical characteristic of a creep crack strain curve is that the slope of the curve decreases with time until stabilizing. Furthermore, the steady-state slope of the creep crack strain curve increases with the loading level. As the loading level increases to a certain value, the rock crack strain rate first decreases, then stabilizes, and ultimately increases due to rock crack development. The axial creep crack strain rate of micritic bioclastic limestone under the triaxial multilevel compression creep test at a confining pressure of 1 MPa is plotted in Fig 15. The creep crack strain rate is large at the beginning of creep and gradually decreases and stabilizes with time. The higher the creep stress level is, the larger the steady creep crack strain rate. Studying the properties of the rock creep crack rate evolution will help to establish the rock creep crack evolution equation from the perspective of the rock failure mechanism. Therefore, it is recommended to focus on the crack growth velocity.
The Charles [46] law for subcritical crack extension in glass at a given stress state presents a power law expression for the microcrack growth velocity. The velocity of crack growth V c is usually found to be uniquely related to the stress intensity factor (for a given set of environmental conditions). Other studies on rocks have found that the crack growth characteristics of rock are similar to those of glass [47]. According to the literature [48], some evidence indicates that there is a threshold below which no crack growth occurs. Considering the threshold effect, the crack growth velocity can be expressed as: To normalize the stress intensity factor difference, Eq (13) can be replaced by the following: where V c denotes the crack growth velocity, A and n denote material parameters, K 0 is the threshold of stress intensity factor for growth, K c denotes the critical value of K i , and <> is the Macaulay bracket. The lack of an accurate initial crack length of the rock makes Eq (14) difficult to use directly because the stress intensity factor is related to the initial crack length. To overcome this difficulty, considering the similarity between microscopic cracks and macroscopic crack strains, the crack growth velocity V c can be expressed by Eq (15) [11,48]. where σ p denotes the peak deviatoric stress and σ ci denotes the crack initial stress, which is the threshold of creep deformation. R d , which characterizes the bearing state of the rock, is introduced to replace the stress in Eq (16) for dimensionless purposes. Considering that the crack growth velocity is proportional to the crack strain rate, the creep crack strain rate can be written as: Further simplification yields the expression for the steady-state creep crack strain rate: Fig 15 shows that the axial creep crack strain rate of the rock gradually decreases with time and eventually stabilizes at the steady creep rate. The growth and adjustment of rock microcracks and pores when the rock is subjected to a constant load result in a gradual decrease in the rock crack growth velocity. The adjustment is expressed macroscopically as a decrease in the rock crack strain rate with time. Therefore, the rock creep crack rate evolution equation based on the steady state creep rate should be written as: where m and B denote two material parameters.
Taking an integral over Eq (18) gives: where F denotes the parameter of loading crack strain and −m×(R ci d −R ci d )/ B+F is the instantaneous loading crack strain.
Previous researchers have proposed elements including springs, sliders, dashpots and fractional order elements that cannot characterise the evolution of rock crack strain over time. It is necessary to develop elements that characterise the evolution of rock creep crack. According to Eq 19, the element representing the creep crack strain of the rock, which is referred to as Mo's element, is proposed, as shown in Fig 16. Based on the creep crack strain evolution element proposed in this paper, the standard particle swarm optimization (SPSO) algorithm was employed to fit the experimental data. The parameters A and n are correlated, i.e., their values affect each other. To solve this problem, simplify the equation and reduce the number of parameters, the parameter A was set to 1. The fitting and test curves are shown in Fig 17, demonstrating that the creep crack strain element can describe the creep strain evolution of the rock well. The model parameters are shown in Table 5. The fitted creep crack strain element parameters n, m, B and F in Table 5

Unified transient creep constitutive model (Mo's model)
Based on the above analysis of the transient and creep crack characteristics of the rock and the creation of crack elements, a unified transient creep constitutive model, which is referred to as In accordance with the crack strain evolution characteristics in triaxial compression tests, combined with R d , the mechanical properties of the rock can be represented by the proposed constitutive model as: 8 > > < > > : Based on Eqs (12), (19) and (20), considering that the instantaneous loading crack strain equals −m×(R ci d −R ci d )/ B+F, the rock transient creep constitutive model can be expressed as: 8 > > > > > > > > > > > > < > > > > > > > > > > > > :

Conclusions
The crack strain evolution characteristics of micritic bioclastic limestone obtained from the No. 1 drainage tunnel of the Artashi Hydropower Project in Xinjiang, China, were studied by performing triaxial multilevel loading creep tests under different confining pressures. The following conclusions were obtained.
1. For the rock volumetric strain curve without a typical reversal arc to determine σ cd , an improved method is proposed based on the axial and lateral crack strain, which can determine the appropriate stress threshold for micritic bioclastic limestone.
2. R d , representing the bearing state of the rock, is introduced to investigate the crack evolution characteristics of the rock in triaxial compression tests. An elastic-crack model is proposed based on R d , and the crack evolution law, test data and model solution fit well with each other, indicating the suitability of this model to describe and predict the transient crack evolution of micritic bioclastic limestone.
3. Although the micritic bioclastic limestone has a high elastic modulus, the high content of bioclast matter makes the creep of the rock significant, which results in long-term deformation and fracture development throughout the life of engineering projects. From creep tests under different confining pressures, a creep crack strain evolution equation was obtained via the integral of the crack evolution rate and fits the experimental data closely.